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Abstract. - We study numerically the influence of density and strain rate on the diffusion and 
mobility of a single tagged particle in a sheared colloidal suspension. We determine independently 
the time-dependent velocity autocorrelation functions and, through a novel method, the response 
functions with respect to a small force. While both the diffusion coefficient and the mobility depend 
on the strain rate the latter exhibits a rather weak dependency. Somewhat surprisingly, we find 
that the initial decay of response and correlation functions coincide, allowing for an interpretation 
in terms of an 'effective temperature'. Such a phenomenological effective temperature recovers the 
Einstein relation in nonequilibrium. We show that our data is well described by two expansions 
to lowest order in the strain rate. 
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Introduction. — The mobility of a single spherical 
particle immersed in a solvent determines the velocity of 
the particle in response to an applied external force. For 
small Reynolds numbers Stokes' law yields the famous ex- 
pression (Iq 1 = 3Trr]a in terms of the sphere diameter a and 
the solvent viscosity r\ in thermal equilibrium. This free 
mobility [1q is intimately related to spontaneous solvent 
fluctuations through the Einstein relation. For a suspen- 
sion of interacting particles, even without hydrodynamic 
coupling, the mobility \i of a single tagged particle is re- 
duced. This reflects the fact that work is necessary to dis- 
place neighboring particles in order for the tagged particle 
to move, leading to larger dissipation. Still, in equilibrium 
the Einstein relation 



D = k^T/i 



(1) 



equates the effective, long-time diffusion coefficient D ob- 
tained from measuring the mean square displacement of a 
single tagged particle with its mobility through the solvent 
temperature T, where fee is the Boltzmann constant. 

The Einstein relation ([!]) is one out of many fluctuation- 
dissipation relations valid in the linear response regime for 
small perturbations of the equilibrium state [lj. It is cru- 
cial to realize that also nonequilibrium steady states allow 
for a linear response. However, driving the suspension 
beyond the linear response regime, fluctuation-dissipation 



relations such as the Einstein relation ([I]) need to be gen- 
eralized to nonequilibrium. There are two basic strategies 
discussed in the literature. The first strategy is to intro- 
duce an additive correction taking on the form of another 
correlation function [2j|4]. This correlation function in- 
volves another observable that can be related either to 
entropy production [5] or 'dynamical activity' j6j. Such 
an approach has been demonstrated experimentally for a 
single driven colloidal particle [7j-[9]. The second strategy 
introduces a multiplicative correction through an effective 
temperature 10 11 replacing T in Eq. ([!]). 



11 replacing T in Eq. (|lj). Originally 
developed in the context of aging, glassy dynamics and 
weakly driven systems, effective temperatures have also 
been investigated in shear driven supercooled liquids [12] . 

Self-diffusion in sheared interacting suspensions has 
been studied extensively in computer simulations 
and experiments [17||18| as well as analytically 
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large body of publications studies supercooled conditions 
in relation with the glass transition 12 18 22 . Most of 



these works focus on the self-diffusion coefficient as this 
quantity is easily obtained from experiments and simu- 
lations. The mobility of a tagged particle has been ad- 
dressed somewhat less prominently and mostly in ana- 



lytic calculations 21 23 . In this Letter we determine nu- 



merically both the full time-dependent velocity response 
and autocorrelation functions for a tagged particle of a 
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Fig. 1: a) Simple shear flow with strain rate 7. b) Pair distribu- 
tion function g(r) in the xy plane for volume fraction (j> — 0.3 
and strain rate 7 = 1. Centered on any particle, the function 
g(r) quantifies the probability to find another particle at r. 
While isotropic in equilibrium, the pair distribution function is 
distorted through external flow. 



hard-core Yukawa suspension driven into a nonequilibrium 
steady state through simple shear flow. In contrast to pre- 
vious work, explicit knowledge of the response function 
allows us to calculate and discuss the mobility as a func- 
tion of density and strain rate. We use a novel method 
to efficiently obtain the time-dependent response function 
of the velocity with respect to a small force applied to a 
single particle. Similar methods to extract the response of 
a system using only unperturbed steady state trajectories 
have been discussed in Refs. |24||25| . 

Sheared hard-core Yukawa fluid. The N col- 
loidal particles interact through the purely repulsive 
Yukawa potential 



u(r) 



00 



(r^l) 
(r < 1) 



(2) 



with hard core exclusion. The two potential parameters 
are the energy e at contact and the screening length re -1 
determining the range of interactions. Changing re in- 
terpolates between hard-sphere (large re) and coulombic 
(small re) interactions. Throughout the paper we employ 
dimensionless units and measure length in units of the 
particle diameter a and energies in units of k^T. The 
time scale 'iira 3 r]/k-QT is set by the time a particle diffuses 
a distance equal to its diameter. In particular, employ- 
ing these units the mobility and diffusion coefficient of 
a free particle reduce to unity, D = fi = 1. We set 



= 0.2 and choose e = 



such that the liquid is sta- 
We explore the liquid 
0.1,0.2,0.3,0.4, where 



20 



ble for a large pressure range 
phase at four volume fractions 
4> = nN(a/L) 3 /6 with L the side length of the cubic sim- 
ulation box. The highest density cj> = 0.4 is close to the 
equilibrium freezing transition, which occurs at a pressure 
28.9 [26] (for <f) = 0.4 the measured pressure in our simu- 
lation is 27.6). For comparison, the freezing transition in 
a hard sphere suspension occurs at <f) ~ 0.494 [27) . 

We employ Brownian dynamics simulations, for de- 
tails see the appendix. The suspension is driven into 
a nonequilibrium steady state through simple shear flow 



u(r) = ^ye x with strain rate 7 (which equals the Peclet 
number in our units), see Fig. (TJl) . The equations of mo- 
tion are r^ = v° and 



-V fc J7 - [v° - u(r fc )] + f fc + £ 



k 1 



(3) 



where the dimensionless mass is set to one. Physically, 
this choice implies that momenta relax on the diffusive 
time scale. Besides the forces due to the potential en- 
ergy U = J2k<k' u (\ T k - r k'\) we allow for direct forces f fc . 
The stochastic noise £ fe modeling the interactions with the 
solvent has zero mean and correlations (£fci(i)£fc'j(i')) = 
2SijSkk'3{t — t'), where i,j = x,y,z is the vector com- 
ponent. In Eq. ([3]), we neglect hydrodynamic coupling 
between different particles due to the solvent. 

For the shear flow turned on we correct the particle 
velocities {v°} by the external flow and investigate the 
relative velocity v/. = v° — u(rfe). We are interested in 
the dynamics of a single tagged particle interacting with 
the remaining TV — 1 particles in the suspension. Since all 
particles are identical we designate particle 1 as the tagged 
particle and drop the subscript; in the following r and v 
are the position and relative velocity of the tagged particle, 
respectively. We define the response of this velocity 



Rij(t - t';j) 



S(vi(t)) 
Sfj(t') 



(4) 



with respect to an additional small force f . In addition, 
we define the relative velocity autocorrelation matrix 



C tJ (t-t';j) = (« i (t)« J -(0)o- 



(5) 



The brackets (• • -)o refer to an average with respect to the 
unperturbed steady state whereas (•) is the average with 
the external force applied. In equilibrium (7 = 0) the 
fluctuation-dissipation theorem Rij(t) — Cy(t) holds. 

Response function. — Sampling the correlation 
function ^ is straightforward. The direct way to ob- 
tain the response matrix Q from simulations would be 
through a step perturbation of the force and the subse- 
quent recording of the tagged particle's velocity. Such a 
protocol has to be repeated many times and the corre- 
sponding response function follows as the time-derivative 
of the mean velocity. Here, we follow another route and 
determine the response function through the path inte- 
gral representation of the fluctuation-dissipation theorem 
(FDT) for nonequilibrium steady states. The FDT in its 
general form reads 

Rijit-t'-n) = (» i (t)B i (t'))o 

with observable Bj(t) — SlnP/Sfj(t)\f—o conjugate to the 
perturbation force f acting on the tagged particle. The 
stochastic path weight is 

P[{t k (t)Y, f(*)] = A^exp j- \ J2 J dt €*(*)} 
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Fig. 2: Comparison of the off-diagonal components R xy (t) and 
R yx (t) for strain rate 7 = 1.0 and the four simulated volume 
fractions: a) <j> = 0.1, b) (/> = 0.2, c) <j> = 0.3, and d) <j> = 0.4. 
Increasing the density, the two curves approach each other until 
for the highest density they almost lie on top of each other. 
Note the changing time scale. 



with normalization constant J\f. From Eq. ([3| we see that 
a perturbation of f is equivalent to a perturbation of £ 
with 6£i(t)/6fj(t') = —SijS(t — t') we immediately obtain 
Bj — £j/2 and therefore 



Rij(t 



(6) 



Since in a computer simulation we have direct access to 
the noise we can exploit such an expression to obtain the 
response function through a steady state correlation func- 
tion. While Eq. ([6| has been known before [4 28 , to the 
best of our knowledge so far it has not been exploited to 
obtain the response function numerically. 

To understand the influence of the shear flow on the 
particle motion it is instructive to look at the off-diagonal 
components R xy and R yx plotted in Fig. [2j The compo- 
nent R yx describes the mean behavior of a tagged particle 
when we apply a force parallel to the shear flow and mea- 
sure its velocity perpendicular in y direction. The behav- 
ior of R yx can be explained by looking at the pair distribu- 
tion function in Fig. [Ij}), which is deformed compared to 
its equilibrium isotropic shape. In order to move faster at 
short times the particle moves up (R yx is positive) to over- 
come its neighbors through a region of lower probability 
to encounter another particle. At a later time the particle 
is pushed back (R yx is negative) due to interactions with 
other particles which become more pronounced at higher 
densities. Interchanging x and y-direction, the same argu- 
ments hold for the component R xy . However, since we pull 
the particle up it enters a region where the surrounding 
fluid moves faster due to the shear flow. Hence, the parti- 
cle is accelerated and the response of the relative velocity 



is negative for small times. With increasing density the 
particle cannot move far in y-direction, making the veloc- 
ity differences smaller. The qualitative difference between 
the two curves diminishes and for 4> — 0.4 both almost lie 
on top of each other. Hence, the effect of the shear flow on 
a single particle is more and more symmetric as particle 
motion becomes correlated at higher densities. 

Diffusion and mobility. We now turn to the 
nonequilibrium diffusion coefficients and mobilities, 



Di 



dtCij(t), fiij 



= 9^) 



AtR l3 {t), (7) 



which are obtained through integrating the velocity au- 
tocorrelation and response matrix, respectively. The mo- 
bility is defined as the velocity change in response to a 
small force applied to the tagged particle, perturbing the 
steady state reached through shearing the solvent. The 
diffusion coefficients are related to the velocity autocor- 
relation through a Green-Kubo kind relation. They are 
plotted in Fig. [3] and increase with increasing strain rate. 
From our data we find that we can distinguish diffusion 
parallel to the shear flow with Dm = D xx and diffusion 
perpendicular with D± — D yy — D zz . In shear flow the 
tagged particle moves between layers of different flow ve- 
locity effectively leading to larger fluctuations, allowing 
the particle to explore phase space faster. At low density 
the diffusion Dm parallel to the shear flow is substantially 
enhanced compared to D± even though we subtract out 
the external flow. However, the difference between the two 
diffusion coefficients vanishes with increasing density. 

In Fig. [4] we plot the reduced mobility //(</>, 7)/ /x e q(0) 
versus strain rate 7 for the four simulated densities with 
equilibrium (7 = 0) mobility /i oq . We find that the diag- 
onal components of the mobility matrix are equal within 
error bars and we obtain the shown mobilities through 
averaging over the three directions. The off-diagonal com- 
ponents are somewhat harder to obtain due to large sta- 
tistical errors but are clearly much smaller than their di- 
agonal counterparts (data not shown). The dependence 




Fig. 3: Diffusion coefficients Du parallel and D± perpendicular 
to the shear flow vs. strain rate 7 for the four different volume 
fractions <f>. For a free particle D±_ = Du = 1. 
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Fig. 4: Reduced mobilities /V/ieq vs - strain rate 7 for the four 
different volume fractions (f>. The lines are fits to Eq. j8J. 



of the absolute value of the mobility on the strain rate is 
rather weak and for the lowest density = 0.1 it is even 
constant. Such a weak dependence suggests that the abil- 
ity of the solvent to reorganize in response to dragging 
the tagged particle out of the 'cage' formed by neighbor- 
ing particles is only slightly affected by the presence of the 
shear flow. Going to supercooled conditions, this is likely 



to break down 29 



To explain the dependence of the mobility on the strain 
rate we consider the mean velocity of the tagged particle 
from Eq. 



f = -{N/Lf J dr j(r)Vn(r) + f . 



Here, <?(r; tf>, 7, f ) is the pair distribution function to find a 
second particle at r if there is a particle at the origin, see 
Fig. [lj>), and = — Vif/ is the force exerted by neigh- 
boring particles on the tagged particle. The effects of the 
shear flow and the force f enter only through the struc- 
ture information encoded in the pair distribution function. 
We can expand g into a Taylor series for small forces f . 
On the other hand, it is well known that the structure of 
the suspension in the presence of shear flow is singularly 

re- 



perturbed from its isotropic equilibrium form 30 31 
quiring an asymptotic expansion in powers of -y 1 ^ 
an expansion in lowest order leads to 



jU((£,7) w MeqOW l + x(</>)7 



Such 



(8) 



In principle the coefficients x{4>) could be obtained from 
the knowledge of the perturbed pair distribution function 
<?(r). Here, we determine them through fitting the mobil- 
ity, see the lines in Fig. [4} The fits show a good agreement 
with the simulated data for all strain rates and densities 
even though we have only retained the lowest order of the 
expansion ([8]). 

Einstein relation. — In Fig. [5] we plot R xx (t) and 
C X x(t) as functions of time for different volume fractions. 
The correlation functions have been scaled by a constant 



factor 1 /6 X to match the initial decay of the response func- 
tions. This procedure reveals a rather good agreement 
between response and correlation function even for longer 
times (this holds also for the yy and zz components). We 
approximate the ratio 



i(0,7) w i + c^O/Ot 2 



(9) 



by these constant factors. Using Ru(0) = 1 we can inter- 
pret 9i fa {v?)o as effective temperatures since they equal 
approximately the velocity fluctuations of the tagged par- 
ticle. We have also checked the distribution functions of 
the velocity which are Gaussian with width \[9i as ex- 
pected. 

In Eq. ([9]) we expand the correlation function to second 
order in the strain rate (due to the symmetry of simple 
shear flow there is no first order) with fit parameters o^. 
This expansion becomes exact in the case of interacting 
particles with linear forces [32] . In the insets of Fig. [5] the 
factors Qi for the three directions x, y, and z are shown to 
follow this quadratic prediction. Similar to the diffusion 
coefficients we can distinguish a factor 6u = 9 X parallel, 
and a factor 9± = 9 y = 9 Z perpendicular to the shear 
flow. Again, increasing the density the difference between 
the two directions vanishes. 

Two points are noteworthy. First, no simple proportion- 
ality can be found between the off-diagonal components of 
the response and the correlation matrix. Both have a qual- 
itatively different shape (data not shown). In particular, 
the off-diagonal response components are strictly zero at 
t = whereas, in nonequilibrium, the off-diagonal veloc- 
ity correlations are different from zero. Second, there is a 
fundamental difference compared to the effective temper- 
ature discussed for non-stationary, glassy dynamics out 
of equilibrium. There, fluctuation and dissipation are re- 
lated by an effective temperature at low frequencies (i.e., 
on long time scales) while the initial decay of response 



and correlations is governed by the bath temperature 10 



Moreover, the effective temperature evolves slowly as the 
system approaches equilibrium. In contrast, in our case 
already the initial decay is governed by a temperature 
6i > 1. The effect of this temperature extends into the 
tails of response and correlation functions. It is only for 
high densities that we observe a deviation in the tails as 
can be seen in Fig. [5j:) . 

We can finally write down a simple generalized Einstein 
relation 



if! 



(10) 



for the diagonal components of the diffusion matrix. In- 
serting the expansions for mobility [Eq. ^} and effective 
temperature [Eq. we find that D^ — /z eq ~ 7 1 / 2 to low- 
est order. Such a dependence was also found in molecular 



dynamics simulations for a Lennard-Jones fluid 13 , 14 
Due to the small x/ a i ratios we cannot resolve this 7^" 
dependence here. 
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Fig. 5: The response R xx (t) (dashed) and correlation C xx (t) (dotted) functions for volume fractions a) <f> = 0.1, b) (j> = 0.2, and 
c) (f> = 0.3. The correlation functions are scaled by a constant factor 1/9 X to match the initial decay of the response function. 
The insets show these factors as a function of strain rate for the three directions, where the lines are fits to Eq. d9). 



In Fig. [6] we test the putative effective temperature by 
comparing 8±/i to the numerically obtained diffusion coef- 
ficient D± perpendicular to the shear flow. The mobility 
for <j) = 0.1 is independent of strain rate and the diffusion 
follows the quadratic prediction D± oc 7 2 . While we find 
a good agreement for the two lowest densities, the effec- 
tive temperature underestimates the diffusion coefficient 
at intermediate strain rates and high densities since the 
diffusion coefficient qualitatively changes and approaches 
a linear function D± oc 7 at high densities. This indicates 
that the differences in the tails of response and velocity 
autocorrelation funtions become more important. Also, 
higher order terms might be required in the expansion of 
mobility and effective temperature. 

Experimental realization. — We briefly discuss how 
our findings could be tested in experiments. Of course, the 
route via Eq. ([6| to obtain the velocity response matrix 
through the explicit knowledge of the noise is not available 
in experiments. Moreover, the direct route, i.e. perturbing 
only a single particle within a suspension and measuring 
its time-dependent mean velocity, is experimentally chal- 
lenging and, as we find in our simulations, also statistically 
more demanding. 




Fig. 6: Test of the Einstein relation D± = 8±n with effective 
temperature #x from Eq. |9| and mobility from Eq. |8j. The 
curves show a very good agreement for the lowest two densities 
but start to deviate at higher densities. 



Despite the difficulties it is still interesting to obtain this 
response function since it immediately yields the nonequi- 
librium mobility. We now assume that the tagged particle 
undergoes overdamped motion which is certainly the rele- 
vant limit for experiments. The Langevin equation for the 
tagged particle reads 



(11) 



The replacement of the noise in Eq. j6) by Eq. (|TTJ) is per- 
missible since the Jacobian arising due to the change of 
variables is independent of f . We then obtain an experi- 
mentally accessible expression for the response function 



(12) 



for components i,j = y,z perpendicular to the shear flow. 
Let us assume that we record the particle position at 
times tk = kr with time resolution r, e.g., through video 
microscopy. The velocity is then approximated through 
the finite difference = (r^ — T^-ij/r. In principle, the 
force on the tagged particle can be calculated from 
the knowledge of the pair potential and the positions of 
all neighboring particles. 

Conclusions. — We have studied a hard-core Yukawa 
colloidal system at different densities driven into a 
nonequilibrium steady state through shear flow. In par- 
ticular, we investigate diffusion and mobility of a single 
tagged particle for four volume fractions <j) and intermedi- 
ate strain rates 7^1. The self-diffusion coefficient is cal- 
culated through the Green-Kubo relation from the single 
particle's velocity autocorrelation function. The mobility 
is obtained from the particle's response function through 
integration. For systems governed by stochastic dynamics, 
this response function can be obtained efficiently from the 
correlation function Eq. ^ measured in the unperturbed 
steady state. While for low densities we can clearly dis- 
tinguish quantities measured parallel and perpendicular to 
the shear flow, this difference vanishes for high densities. 

Surprisingly, the diagonal components of the response 
(i.e., the response is measured in the direction of the ap- 
plied force) and correlation matrix can be matched over 
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a large time range. The resulting proportionality factor 
can be interpreted as an effective temperature, effectively 
restoring the Einstein relation connecting diffusion and 
mobility. Moreover, this proportionality factor is well ap- 
proximated by a quadratic expansion in the strain rate. 
It will be important to study how general such a simple 
effective temperature is for driven interacting colloidal sus- 
pensions and whether it extends to other observables. We 
believe that the methodology presented here will lead to 
new insights in the numerical study of dense colloidal sus- 
pension, e.g., for microscopic stress fluctuations 32 35 



Finally, the influence of hydrodynamic interactions on our 
results remains to be investigated. 
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Appendix: Simulation details. — The simulated 
systems consist of N = 1728 particles in a cubic simu- 
lation box. Since we are interested in the bulk behavior 
of the suspension we employ periodic boundary conditions 
and account for the shear flow through Lees-Edwards slid- 
ing bricks. The equations of motion are integrated by 
a stochastic version of the velocity Verlet algorithm [33] , 
where the velocity appearing in the force term at the right 
hand side of Eq. ([3| is taken from the mid-step velocity. 
The time step is set to At = 0.0005 < (kv^)' 1 ~ 
0.08 ... 0.2. We equilibrate the suspension and then slowly 
increase the strain rate to the final value. Correlation func- 
tions have been obtained by averaging over 400 particle 
trajectories with 300,000 time steps each. These trajecto- 
ries have been determined in two independent runs from 
randomly chosen particles. 

To implement the hard-core repulsion and prevent par- 
ticles from overlapping, the following simple algorithm is 
employed (see also Refs. 16 34 and references therein). 



After every particle has been moved, but before new forces 
are calculated, we store all overlapping particle pairs and 
remove these overlaps as follows. For each pair both of the 
particles are moved backwards in time along their respec- 
tive velocity vector up to the point where they collided. 
This time < s < At is stored. Knowing the positions 
and velocities at the impact, we compute the connection 
vector e between both particles. We decompose the ve- 



locities into v 



T 

ee J vi 2 



and 



ee T )vi 2 - 



'1,2 — cc y i^ mlu v i,2 — \ 
Only the parts parallel to e can change. Using the usual 
elastic collision rule preserving momentum and kinetic en- 
ergy of the particles we obtain the after-collision velocities 



Bechinger, Phys. 
Phys. Rev. 



and Vn = 



From the positions of 



their collision the particles are then propagated forward 
with time step s along the new velocity vectors. This pro- 
cedure is repeated as long as overlapping pairs exist. 
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